rm(list = ls())
europeLink="https://github.com/Spatial-Data-Analytics-DACSS-690D/HW1-Interactive-Visualization/raw/refs/heads/main/WORLD/europe_8857.gpkg"
library(sf)
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
st_layers(europeLink)
europe <- st_read(europeLink, layer = "continent")
## Reading layer `continent' from data source
## `https://github.com/Spatial-Data-Analytics-DACSS-690D/HW1-Interactive-Visualization/raw/refs/heads/main/WORLD/europe_8857.gpkg'
## using driver `GPKG'
## Simple feature collection with 37 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -2668477 ymin: 3480747 xmax: 3207182 ymax: 7761460
## Projected CRS: WGS 84 / Equal Earth Greenwich
europemob_ddm <- st_read(europeLink, layer = "mobiles_ddm")
## Reading layer `mobiles_ddm' from data source
## `https://github.com/Spatial-Data-Analytics-DACSS-690D/HW1-Interactive-Visualization/raw/refs/heads/main/WORLD/europe_8857.gpkg'
## using driver `GPKG'
## Simple feature collection with 14480 features and 0 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -2662224 ymin: 3515365 xmax: 3197662 ymax: 7749696
## Projected CRS: WGS 84 / Equal Earth Greenwich
europemob_psm <- st_read(europeLink, layer = "mobiles_psm")
## Reading layer `mobiles_psm' from data source
## `https://github.com/Spatial-Data-Analytics-DACSS-690D/HW1-Interactive-Visualization/raw/refs/heads/main/WORLD/europe_8857.gpkg'
## using driver `GPKG'
## Simple feature collection with 37 features and 5 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -1290663 ymin: 4371225 xmax: 2918367 ymax: 7326054
## Projected CRS: WGS 84 / Equal Earth Greenwich
st_crs(europe)$epsg
## [1] 8857
library(ggplot2)
library(scales)
ggplot() + theme_light() +
#call map
geom_sf(data = europe) +
# adjust datum
coord_sf(datum = st_crs(europe)) +
# custom X and Y axes
scale_x_continuous(labels = label_number(scale = 1/1000, suffix = " km")) +
scale_y_continuous(labels = label_number(scale = 1/1000, suffix = " km")) +
labs(x = "Easting", y = "Northing") #labels
# DDM - static
plot1_a=ggplot() + theme_light() +
#base
geom_sf(data = europe,fill="white",color="grey80") +
# layer
geom_sf(data=europemob_ddm,
alpha=0.1,#transparency
shape=20,
size=0.5,
color="red")+
coord_sf(datum = st_crs(europe)) +
# custom X and Y axes
scale_x_continuous(labels = label_number(scale = 1/1000, suffix = " km")) +
scale_y_continuous(labels = label_number(scale = 1/1000, suffix = " km")) +
labs(x = "Easting", y = "Northing") #labels
plot1_a
# DDM - interactive
library(leaflet)
# needs geographic CRS
ddm_4326 <- st_transform(europemob_ddm, crs = 4326)
leaflet() %>%
addTiles() %>% # base map
addCircleMarkers(data = ddm_4326,
color="red",
stroke=0,
radius = 0.01,
opacity = 0.1)
Almost the same:
plot1_b=leaflet() %>%
addProviderTiles(providers$CartoDB.Positron) %>% # names in english
addCircleMarkers(data = ddm_4326,
color="red",
stroke=0,
radius = 0.01,
opacity = 0.1)
plot1_b
plot2_a=ggplot() + theme_light() +
#base
geom_sf(data = europe,fill="white",color="grey80") +
# layer
geom_sf(data=europemob_psm,
shape=20,
aes(size=mobiles/1000,color=factor(mobiles_outlier)))+
scale_size_area(
name = "Mobile Users (x1,000)", # Name for the legend
max_size = 20) + # This defines the maximum *radius*
guides(color="none") +
coord_sf(datum = st_crs(europe)) +
# custom X and Y axes
scale_x_continuous(labels = label_number(scale = 1/1000, suffix = " km")) +
scale_y_continuous(labels = label_number(scale = 1/1000, suffix = " km")) +
labs(x = "Easting", y = "Northing") #labels
plot2_a
# PSM - interactive
# CRS for leaflet
psm_4326 <- st_transform(europemob_psm, crs = 4326)
plot2_b=leaflet() %>%
addProviderTiles(providers$CartoDB.Positron) %>%
addCircleMarkers(data = psm_4326,
# fillColor = "red",
# stroke = T, #border
# color="red", #border color
# fillOpacity = 1,
label = ~Country,
radius = ~size)
plot2_b
plot3_a=ggplot() + theme_light() +
#base
geom_sf(data=europe,
aes(fill=mobiles_density_FJ5_cat))+
coord_sf(datum = st_crs(europe)) +
scale_fill_brewer(direction = 1,palette = 'RdBu') +
scale_x_continuous(labels = label_number(scale = 1/1000, suffix = " km")) +
scale_y_continuous(labels = label_number(scale = 1/1000, suffix = " km")) +
labs(x = "Easting", y = "Northing",fill="Mobile Users Density")
plot3_a
choro_4326 <- st_transform(europe, crs = 4326)
#needed
choro_4326$mobiles_density_FJ5_cat=ordered(choro_4326$mobiles_density_FJ5_cat)
# First, create a color palette for your categorical variable
pal <- colorFactor(palette = "RdBu",
domain = choro_4326$mobiles_density_FJ5_cat,
ordered = T,
reverse = F)
leaflet() %>%
addProviderTiles(providers$CartoDB.Positron) %>%
addPolygons(
data = choro_4326,
fillColor = ~pal(mobiles_density_FJ5_cat),
fillOpacity = 0.7,
color = "white",
weight = 1,
opacity = 1
) %>%
addLegend(
pal = pal,
values = choro_4326$mobiles_density_FJ5_cat,
)
# Calculate the center of europe for the home button
europe_center <- st_bbox(choro_4326) %>%
as.vector() %>%
{c(mean(c(.[1], .[3])), mean(c(.[2], .[4])))} # Calculate center coordinates
plot3_b=leaflet() %>%
addProviderTiles(providers$CartoDB.Positron) %>%
addPolygons(
data = choro_4326,
fillColor = ~pal(mobiles_density_FJ5_cat),
fillOpacity = 0.7,
color = "white",
weight = 1,
opacity = 1,
label = lapply(
paste("<strong>", choro_4326$Country, "</strong><br>",
"Density: ", choro_4326$mobiles_density),
htmltools::HTML
)
) %>%
addLegend(
pal = pal,
values = choro_4326$mobiles_density_FJ5_cat,
opacity = 0.7,
title = "Mobile Users Density",
position = "bottomright"
)%>%
addEasyButton(
easyButton(
icon = "fa-home",
title = "Reset to europe View",
onClick = JS(sprintf("function(btn, map){ map.setView([%f, %f], 3); }",
europe_center[2], europe_center[1]))
)
)
plot3_b
saveRDS(plot1_a,file = "plot1_a.rds")
saveRDS(plot1_b,file = "plot1_b.rds")
saveRDS(plot2_a,file = "plot2_a.rds")
saveRDS(plot2_b,file = "plot2_b.rds")
saveRDS(plot3_a,file = "plot3_a.rds")
saveRDS(plot3_b,file = "plot3_b.rds")